Ejercicio de Machine Learning: Clasificación y regresión vinícola

Josep López Lizarte - Máster Data Science

En este trabajo trataremos un ejercicio de machine learning con el objetivo de generar, entrenar, validar y testear modelos tanto de clasificación como de regresión.

Se utilizará un dataset que proviene de la Universidad de Minho, generado por P. Cortez. Dicho dataset se encuentra en el UC Irvine Machine Learning Repository (aquí) está disponible.

A partir del mismo, centraremos nuestro análisis de machine learning:

En primer lugar, importaremos nuestro dataset para iniciar un análisis descriptivo del mismo.

In [116]:
# Importamos bibliotecas necesarias:
import numpy as np
import pandas as pd
import os
%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import pandas_profiling
from IPython.display import Image

Revisamos que el directorio de trabajo del cual queremos extraer el dataset es el correcto, por si fuera necesario modificarlo.

In [2]:
print("Directorio de trabajo actual ", os.getcwd())
Directorio de trabajo actual  C:\Users\34662\Escritorio\MÁSTER DATA SCIENCE\Material Máster Data Science\3 - Python\Prácticas\Practica Ejercicio Machine Learning (clasificacion y regresion)

Con la ayuda de la función read_csv(), tendremos nuestros datos cargados en un dataframe de pandas.

In [3]:
df_vino = pd.read_csv(r'C:\Users\34662\Escritorio\MÁSTER DATA SCIENCE\Material Máster Data Science\3 - Python\Prácticas\Practica Ejercicio Machine Learning (clasificacion y regresion)\winequality.csv',
                      sep = ';')

df_vino.head()
Out[3]:
fixed acidity volatile acidity citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
0 5.20 0.34 0.00 1.8 0.050 27.0 63.0 0.99160 3.68 0.79 14.0 6 red
1 6.20 0.55 0.45 12.0 0.049 27.0 186.0 0.99740 3.17 0.50 9.3 6 white
2 7.15 0.17 0.24 9.6 0.119 56.0 178.0 0.99578 3.15 0.44 10.2 6 white
3 6.70 0.64 0.23 2.1 0.080 11.0 119.0 0.99538 3.36 0.70 10.9 5 red
4 7.60 0.23 0.34 1.6 0.043 24.0 129.0 0.99305 3.12 0.70 10.4 5 white
In [4]:
df_vino["color"].value_counts()
Out[4]:
white    4898
red      1599
Name: color, dtype: int64

El dataset contiene la información sobre distintos vinos, resumiendo las características de los mismos para poder extraer posteriormente si se trata de un vino tinto (rojo) o blanco. El número de entradas es de 1599 para el vino tinto y 4898 para el vino blanco.

O bien para extraer a través de regresiones la calidad del vino con valoraciones de 0 a 10.

A nivel de estructura, los atributos siguientes representan los inputs del dataset, ya que son los que están basados en tests fisicoquímicos por expertos que han analizado un muestreo de vinos:

1 - fixed acidity

2 - volatile acidity

3 - citric acid

4 - residual sugar

5 - chlorides

6 - free sulfur dioxide

7 - total sulfur dioxide

8 - density

9 - pH

10 - sulphates

11 - alcohol

En cambio, la variable output es la calidad (quality) del vino, que tiene una puntuación de entre 0 y 10 dependiendo de la calidad, y su variación se debe a la sensibilidad sobre las variables descritas anteriormente.

La variable de color será la que tendremos que ser capaz de clasificar lo mejor posible a lo largo del trabajo, por lo que a partir del resto de variables debemos ser capaces de clasificar el color de los vinos.

Posteriormente procederemos a realizar un modelo regresor que prediga lo mejor posible la calidad de los vinos, tratandola como el output de nuestro dataset.

In [5]:
df_vino.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 6497 entries, 0 to 6496
Data columns (total 13 columns):
 #   Column                Non-Null Count  Dtype  
---  ------                --------------  -----  
 0   fixed acidity         6497 non-null   float64
 1   volatile acidity      6497 non-null   float64
 2   citric acid           6497 non-null   float64
 3   residual sugar        6497 non-null   float64
 4   chlorides             6497 non-null   float64
 5   free sulfur dioxide   6497 non-null   float64
 6   total sulfur dioxide  6497 non-null   float64
 7   density               6497 non-null   float64
 8   pH                    6497 non-null   float64
 9   sulphates             6497 non-null   float64
 10  alcohol               6497 non-null   float64
 11  quality               6497 non-null   int64  
 12  color                 6497 non-null   object 
dtypes: float64(11), int64(1), object(1)
memory usage: 660.0+ KB

Con la función .info() vemos que disponemos de un dataset con, aproximadamente 6497 elementos donde cada entrada es un vino. Se aprecia que no es un dataset de un tamaño muy elevado.

Estadística descriptiva

Para poder trabajar con el dataset, es importante conocer que tipo de datos que lo componen, su estructura y forma. Y sobretodo, si existen correlaciones entre variables con las que posteriormente podríamos jugar para obtener la mejor predicción de nuestro modelo.

Realizamos el método .describe() para conocer los valores y estadísticos principales de nuestros atributos.

In [6]:
df_vino.describe()
Out[6]:
fixed acidity volatile acidity citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality
count 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000 6497.000000
mean 7.215307 0.339666 0.318633 5.443235 0.056034 30.525319 115.744574 0.994697 3.218501 0.531268 10.491801 5.818378
std 1.296434 0.164636 0.145318 4.757804 0.035034 17.749400 56.521855 0.002999 0.160787 0.148806 1.192712 0.873255
min 3.800000 0.080000 0.000000 0.600000 0.009000 1.000000 6.000000 0.987110 2.720000 0.220000 8.000000 3.000000
25% 6.400000 0.230000 0.250000 1.800000 0.038000 17.000000 77.000000 0.992340 3.110000 0.430000 9.500000 5.000000
50% 7.000000 0.290000 0.310000 3.000000 0.047000 29.000000 118.000000 0.994890 3.210000 0.510000 10.300000 6.000000
75% 7.700000 0.400000 0.390000 8.100000 0.065000 41.000000 156.000000 0.996990 3.320000 0.600000 11.300000 6.000000
max 15.900000 1.580000 1.660000 65.800000 0.611000 289.000000 440.000000 1.038980 4.010000 2.000000 14.900000 9.000000

Con esto únicamente conoceremos el tipo de números y estadísticos de nuestro dataset, siempre es útil para describir, conocer los valores con los que trabajamos por encima.

Podemos conocer a través del resumen estadístico de la distribución de nuestras features numéricas:

  • Count: el número de valores que no son NaN, en nuestro caso queda reflejado que no hay ninguno, por lo que no será necesario tratar ese apartado.
  • Mean: calcula la media de los valores para cada feature (se ignoran los Nan).
  • Min: nos indica el valor mínimo para cada feature.
  • Max: nos indica el valor máximo para cada feature.
  • Std: muestra la desviación estándar proporcionando una medida de la dispersión de los valores con respecto a la media (mean).
  • 25%, 50%, 75: los cuartiles representan el valor por debajo del cual se situan un determinado porcentaje de los valores para una determinada feature.

La librería de panda profiling lleva varios niveles más allá la información dada por el método describe() de pandas. Generando un informe en HTML y CSS3 que aporta información estadística detallada y gráfica para cada variable de nuestro dataset.

In [7]:
# Guardamos el informe en nuestra carpeta del trabajo.
pR = pandas_profiling.ProfileReport(df_vino)
pR.to_file(output_file="Informe.html")

pandas_profiling.ProfileReport(df_vino)







Out[7]:

Ahora disponemos de un informe analítico con mayores indicativos sobre nuestras variables. Se archiva en la carpeta de trabajo, en formato .html, para ojearlo siempre que lo deseemos

El cual nos realiza una descripción global, nos indica los estadísticos de cada una de las variables, las correlaciones e interacciones que hay entre cada par de variables, si existen datos perdidos e incluso, nos informa de si existen duplicados y los más frecuentes (más adelante trataremos los duplicados).

Es conveniente tener un gráfico de mayor tamaño para poder analizar las variables, por lo que realizamos un histograma de todas las variables para ver la distribución de los datos.

In [8]:
print("Matplotlib version: " + matplotlib.__version__)
Matplotlib version: 3.2.2
In [9]:
df_vino.hist(bins=50, figsize=(20,15))
plt.show()

De este modo, concluímos lo siguiente acerca de las distribuciones de nuestras variables:

1) Nuestras distribuciones principalmente son asimétricas hacia la derecha, por lo que la asimetría es negativa. La media caerá a la derecha de nuestra mediana, por lo que si existe asimetría puede dificultarse para nuestre modelo el hecho de localizar patrones en comparación a como lo haría con una distribución simétrica.

2) No parece que existan outliers relevantes o visualmente destacables. Por lo que no eliminaremos datos de nuestro dataset por el momento para nuestros modelos.

De igual manera podemos realizar un pairplot de Seaborn, así se puede apreciar visualmente con los scatterplots las correlaciones entre ciertas variables. De todas formas al tratarse de un dataset con tantas variables, a nivel visual no se podrá apreciar correctamente, por lo que no podremos extraer conclusiones precisas con el gráfico.

En este caso los colores más rojos son el vino blanco (white) y los mas azules el vino tinto (red).

In [10]:
sns.pairplot(df_vino,
             hue= ("color"),
             vars=[("fixed acidity"),
                   ("volatile acidity"),
                   ("citric acid"), 
                   ("residual sugar"),
                   ("chlorides"),
                   ("free sulfur dioxide"),
                   ("total sulfur dioxide"),
                   ("density"),
                   ("pH"),
                   ("sulphates"),
                   ("alcohol")],
             diag_kind="hist",
             palette="coolwarm",
             height=2
            )
pass

Como se puede apreciar gráficamente, no parece que haya fuertes correlaciones entre las variables, de todas formas verlo gráficamente cuando disponemes de tantas variables es algo complicado a simple vista.

Habrá que verlo mejor con la función .corr() que nos indice el coeficiente de correlación lineal de cada uno de los pares de variables, para cotejar si existe o no correlación entre alguna de ellas.

In [11]:
df_vino.corr()
Out[11]:
fixed acidity volatile acidity citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality
fixed acidity 1.000000 0.219008 0.324436 -0.111981 0.298195 -0.282735 -0.329054 0.458910 -0.252700 0.299568 -0.095452 -0.076743
volatile acidity 0.219008 1.000000 -0.377981 -0.196011 0.377124 -0.352557 -0.414476 0.271296 0.261454 0.225984 -0.037640 -0.265699
citric acid 0.324436 -0.377981 1.000000 0.142451 0.038998 0.133126 0.195242 0.096154 -0.329808 0.056197 -0.010493 0.085532
residual sugar -0.111981 -0.196011 0.142451 1.000000 -0.128940 0.402871 0.495482 0.552517 -0.267320 -0.185927 -0.359415 -0.036980
chlorides 0.298195 0.377124 0.038998 -0.128940 1.000000 -0.195045 -0.279630 0.362615 0.044708 0.395593 -0.256916 -0.200666
free sulfur dioxide -0.282735 -0.352557 0.133126 0.402871 -0.195045 1.000000 0.720934 0.025717 -0.145854 -0.188457 -0.179838 0.055463
total sulfur dioxide -0.329054 -0.414476 0.195242 0.495482 -0.279630 0.720934 1.000000 0.032395 -0.238413 -0.275727 -0.265740 -0.041385
density 0.458910 0.271296 0.096154 0.552517 0.362615 0.025717 0.032395 1.000000 0.011686 0.259478 -0.686745 -0.305858
pH -0.252700 0.261454 -0.329808 -0.267320 0.044708 -0.145854 -0.238413 0.011686 1.000000 0.192123 0.121248 0.019506
sulphates 0.299568 0.225984 0.056197 -0.185927 0.395593 -0.188457 -0.275727 0.259478 0.192123 1.000000 -0.003029 0.038485
alcohol -0.095452 -0.037640 -0.010493 -0.359415 -0.256916 -0.179838 -0.265740 -0.686745 0.121248 -0.003029 1.000000 0.444319
quality -0.076743 -0.265699 0.085532 -0.036980 -0.200666 0.055463 -0.041385 -0.305858 0.019506 0.038485 0.444319 1.000000

Fijandonos bien, no existe ninguna correlación entre ningún par de variables de nuestro dataset que dispongan de una correlación muy fuerte.

Medimos la correlación con el coeficiente de correlación lineal (R), el cual nos lo facilita la función .corr() e interpreta los valores de la covarianza. Nos proporciona valores de entre -1 y 1.

Podemos considerar cierta correlación con las variables siguientes:

  • Density y Alcohol. Donde la correlación nos da un valor de 0.720934.
  • Free sulfur dioxide y Total sulfur dioxide. Siendo una correlación de -0.686745.
  • Density y Residual sugar. Siendo una correlación de 0.552517, muy baja para considerarla significativa.

Las demas variables tienen coeficientes de correlación lineal por debajo de 0.5, por lo que no seria considerable tenerlas en cuenta como variables correlacionadas entre si, es decir, que lo que una varie explique la variabilidad de la otra en esa proporción.

Para cerrar el análisis descriptivo de nuestro dataset, podemos realizar un heatmap de Seaborn.

In [12]:
plt.figure(figsize=(15,10))
sns.heatmap(df_vino.corr(),
            vmin=-1.0,
            vmax=1.0,
            annot=True,
            cmap="coolwarm")
pass

De esta manera se puede contemplar la correlación lineal entre variables de una manera más senzilla y práctica, donde se reflejan con el relleno del número en blanco las correlaciones entre pares de variables más signifcativas (las mismas indicadas más arriba).

Llegamos a la conclusión de que para nuestro ejercicio lo mejor será no eliminar variables, ya que por lo general la mayoría se mantienen en un rango de correlación muy bajo y no sabemos cuanto podría afectar a nuestro modelo negativamente el eliminar variables. Es conveniente realizar primero un modelo predictivo con todas la variables que disponemos, ya que no es un número elevado y no tendremos problema alguno.

Limpieza del dataset

Es importante analizar en primer lugar nuestro dataset, para conocer mejor con que tipo de datos estamos trabajando, pero, posteriormente, es igual de importante revisar los datos y cotejar que esté en las mejores condiciones para tratarlo y poder realizar un modelo de regresión y clasificación con el mismo.

Procedemos a revisar si existen duplicados de lineas en nuestro dataset, para ello utilizamos la función .duplicated().

In [13]:
df_vino[df_vino.duplicated()]
Out[13]:
fixed acidity volatile acidity citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
264 6.6 0.200 0.27 10.9 0.038 29.0 130.0 0.99496 3.11 0.44 10.5 7 white
324 5.2 0.155 0.33 1.6 0.028 13.0 59.0 0.98975 3.30 0.84 11.9 8 white
382 7.7 0.390 0.28 4.9 0.035 36.0 109.0 0.99180 3.19 0.58 12.2 7 white
392 7.2 0.310 0.35 7.2 0.046 45.0 178.0 0.99550 3.14 0.53 9.7 5 white
415 6.8 0.220 0.35 17.5 0.039 38.0 153.0 0.99940 3.24 0.42 9.0 6 white
... ... ... ... ... ... ... ... ... ... ... ... ... ...
6486 6.7 0.410 0.34 9.2 0.049 29.0 150.0 0.99680 3.22 0.51 9.1 5 white
6488 6.9 0.150 0.28 4.4 0.029 14.0 107.0 0.99347 3.24 0.46 10.4 8 white
6489 6.0 0.290 0.25 1.4 0.033 30.0 114.0 0.98794 3.08 0.43 13.2 6 white
6491 6.8 0.270 0.22 17.8 0.034 16.0 116.0 0.99890 3.07 0.53 9.2 5 white
6496 7.1 0.340 0.28 2.0 0.082 31.0 68.0 0.99694 3.45 0.48 9.4 5 red

1177 rows × 13 columns

Una vez detectadas un total de 1177 lineas que están duplicadas, las eliminaremos de nuestro dataset ya que no nos aportan mayor valor a los modelos, son muestras idénticas.

In [14]:
df_vino1 = df_vino.drop_duplicates()
In [15]:
df_vino1[df_vino1.duplicated()]
Out[15]:
fixed acidity volatile acidity citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color

Hemos comprobado que en nuestro nuevo dataset (df_vino1) ya no existen duplicados, por lo que podemos continuar con la limpieza del dataset.

In [16]:
df_vino1.info()
<class 'pandas.core.frame.DataFrame'>
Int64Index: 5320 entries, 0 to 6495
Data columns (total 13 columns):
 #   Column                Non-Null Count  Dtype  
---  ------                --------------  -----  
 0   fixed acidity         5320 non-null   float64
 1   volatile acidity      5320 non-null   float64
 2   citric acid           5320 non-null   float64
 3   residual sugar        5320 non-null   float64
 4   chlorides             5320 non-null   float64
 5   free sulfur dioxide   5320 non-null   float64
 6   total sulfur dioxide  5320 non-null   float64
 7   density               5320 non-null   float64
 8   pH                    5320 non-null   float64
 9   sulphates             5320 non-null   float64
 10  alcohol               5320 non-null   float64
 11  quality               5320 non-null   int64  
 12  color                 5320 non-null   object 
dtypes: float64(11), int64(1), object(1)
memory usage: 581.9+ KB

Para evitar que haya huecos en el índice de nuestro dataset, que no aparezcan saltos de número al haber eliminado los duplicados, resetearemos el índice.

In [17]:
df_vino1.reset_index(drop=True, inplace=True)

Ahora disponemos de un dataset sin duplicados ni huecos en nuestro índice.

Una vez limpiado el dataset, vamos a ordenarlo jerárquicamente para dividir lo que son las features de los labels, así, más adelante será todo más senzillo de codificar y ganaremos tiempo:

In [18]:
df_vino2 = pd.MultiIndex.from_tuples([("features", "fixed acidity"),
                                          ("features", "volatile acitidy"),
                                          ("features", "citric acid"),
                                          ("features", "residual sugar"),
                                          ("features", "chlorides"), 
                                          ("features", "free sulfur dioxide"),
                                          ("features", "total sulfur dioxide"),
                                          ("features", "density"),
                                          ("features", "pH"),
                                          ("features", "sulphates"),
                                          ("features", "alcohol"),
                                          ("label1", "quality"),
                                          ("label2", "color")]
                                        )
df_vino1.columns = df_vino2

Comprobamos nuestra nueva clasificación jerárquica separando entre features y labels:

In [19]:
df_vino1.head()
Out[19]:
features label1 label2
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
0 5.20 0.34 0.00 1.8 0.050 27.0 63.0 0.99160 3.68 0.79 14.0 6 red
1 6.20 0.55 0.45 12.0 0.049 27.0 186.0 0.99740 3.17 0.50 9.3 6 white
2 7.15 0.17 0.24 9.6 0.119 56.0 178.0 0.99578 3.15 0.44 10.2 6 white
3 6.70 0.64 0.23 2.1 0.080 11.0 119.0 0.99538 3.36 0.70 10.9 5 red
4 7.60 0.23 0.34 1.6 0.043 24.0 129.0 0.99305 3.12 0.70 10.4 5 white

Para no perder nuestro dataset en caso de que pueda pasar algo con el ordenador o problemas externos, lo guardaremos a través de la función "pickle" en nuestro disco duro como fichero binario. De este modo nos aseguraremos que no perderemos nunca el dataset y lo podremos cargar cuando lo necesitemos.

In [20]:
df_vino1.to_pickle("vinoteca.pickle")
In [21]:
# Eliminamos el dataset para comprobar que funciona el pickle.
del df_vino1
In [22]:
# Comprobamos que no existe nuestro dataset
try:
    df_vino1.head()
except: NameError

Indicamos que salte el error en caso de producirse para poder ejecutar todas las celdas de código de nuestro Notebook sin problemas.

In [23]:
# Volvemos a generar el dataset a traves de la llamada a nuestro disco duro.
df_vino1 = pd.read_pickle("vinoteca.pickle")
In [24]:
# Vemos que todo está tal como lo dejamos
df_vino1.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5320 entries, 0 to 5319
Data columns (total 13 columns):
 #   Column                            Non-Null Count  Dtype  
---  ------                            --------------  -----  
 0   (features, fixed acidity)         5320 non-null   float64
 1   (features, volatile acitidy)      5320 non-null   float64
 2   (features, citric acid)           5320 non-null   float64
 3   (features, residual sugar)        5320 non-null   float64
 4   (features, chlorides)             5320 non-null   float64
 5   (features, free sulfur dioxide)   5320 non-null   float64
 6   (features, total sulfur dioxide)  5320 non-null   float64
 7   (features, density)               5320 non-null   float64
 8   (features, pH)                    5320 non-null   float64
 9   (features, sulphates)             5320 non-null   float64
 10  (features, alcohol)               5320 non-null   float64
 11  (label1, quality)                 5320 non-null   int64  
 12  (label2, color)                   5320 non-null   object 
dtypes: float64(11), int64(1), object(1)
memory usage: 540.4+ KB

Como podemos apreciar nuestro dataset no dispone de valores NaN (nulos), por lo que no es necesario tratarlos.

Observamos la distribución de nuestra variable categórica, donde la enorme mayoríra son vinos blancos.

In [25]:
df_vino1["label2","color"].value_counts()
Out[25]:
white    3961
red      1359
Name: (label2, color), dtype: int64

Separación de datos en conjunto de train y test

Por el momento tenemos que resolver un problema de clasificación, para saber de la mejor manera posible si un vino es blanco o rojo a partir del resto de variables.

Lo realizaremos con nuestro dataset limpiado de duplicados, pero para el cual no hemos descartado ninguna variable ni hemos realizado reducción de dimensionalidad dado que consideramos por ahora que no es un dataset suficientemente elevado como para hacernos las cosas complicadas.

Es el momento de separar nuestro conjunto de datos en un conjunto de test y otro de training. Nuestro conjunto de test no se utilizará ni tocará hasta que toque evaluar el error de generalización del modelo elegido. Esto nos indica que nuestro trabajo ahora será entrenar el conjunto de training para escoger el mejor modelo de predicción y así, posteriormente, aplicarlo a nuestro conjunto de test para conocer la precisión de acierto del modelo seleccionado.

Dado que sklearn no entiende nada que no sean números, hay que convertir la columna "color" en valores numéricos. Por lo que el 0 será el vino tinto (red) y el 1 el vino blanco (white).

In [26]:
df_vino1.head()
Out[26]:
features label1 label2
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
0 5.20 0.34 0.00 1.8 0.050 27.0 63.0 0.99160 3.68 0.79 14.0 6 red
1 6.20 0.55 0.45 12.0 0.049 27.0 186.0 0.99740 3.17 0.50 9.3 6 white
2 7.15 0.17 0.24 9.6 0.119 56.0 178.0 0.99578 3.15 0.44 10.2 6 white
3 6.70 0.64 0.23 2.1 0.080 11.0 119.0 0.99538 3.36 0.70 10.9 5 red
4 7.60 0.23 0.34 1.6 0.043 24.0 129.0 0.99305 3.12 0.70 10.4 5 white
In [27]:
from sklearn.preprocessing import LabelEncoder

# Creamos instancia
label_encoder = LabelEncoder()

# Con el método .fit() identifica las clases que hay (solo hay 2)
label_encoder.fit(df_vino1["label2","color"])

# Transformamos nuestra columna
df_vino1["label2","color"] = label_encoder.transform(df_vino1["label2","color"])

df_vino1.head()
Out[27]:
features label1 label2
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
0 5.20 0.34 0.00 1.8 0.050 27.0 63.0 0.99160 3.68 0.79 14.0 6 0
1 6.20 0.55 0.45 12.0 0.049 27.0 186.0 0.99740 3.17 0.50 9.3 6 1
2 7.15 0.17 0.24 9.6 0.119 56.0 178.0 0.99578 3.15 0.44 10.2 6 1
3 6.70 0.64 0.23 2.1 0.080 11.0 119.0 0.99538 3.36 0.70 10.9 5 0
4 7.60 0.23 0.34 1.6 0.043 24.0 129.0 0.99305 3.12 0.70 10.4 5 1

Ahora todas las variables de nuestro dataset tienen datos numéricos. Comprobamos que lo ha hecho:

In [28]:
df_vino1.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5320 entries, 0 to 5319
Data columns (total 13 columns):
 #   Column                            Non-Null Count  Dtype  
---  ------                            --------------  -----  
 0   (features, fixed acidity)         5320 non-null   float64
 1   (features, volatile acitidy)      5320 non-null   float64
 2   (features, citric acid)           5320 non-null   float64
 3   (features, residual sugar)        5320 non-null   float64
 4   (features, chlorides)             5320 non-null   float64
 5   (features, free sulfur dioxide)   5320 non-null   float64
 6   (features, total sulfur dioxide)  5320 non-null   float64
 7   (features, density)               5320 non-null   float64
 8   (features, pH)                    5320 non-null   float64
 9   (features, sulphates)             5320 non-null   float64
 10  (features, alcohol)               5320 non-null   float64
 11  (label1, quality)                 5320 non-null   int64  
 12  (label2, color)                   5320 non-null   int32  
dtypes: float64(11), int32(1), int64(1)
memory usage: 519.7 KB

Una vez que todo nuestro dataframe es numérico ya podemos separar en train y test (el 75% del dataset se utilizará para training y el 25% para testing):

In [29]:
from sklearn.model_selection import train_test_split

train, test = train_test_split(df_vino1,
                               train_size = 0.75,
                               test_size = 0.25,
                               random_state = 24)   # Indicamos una seed para no tener variaciones

Observamos que se han generado:

In [30]:
train[:10]
Out[30]:
features label1 label2
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
4261 7.1 0.59 0.00 2.2 0.078 26.0 44.0 0.99522 3.42 0.68 10.8 6 0
5068 7.3 0.25 0.28 1.5 0.043 19.0 113.0 0.99338 3.38 0.56 10.1 6 1
4248 7.4 0.20 0.43 7.8 0.045 27.0 153.0 0.99640 3.19 0.55 9.0 7 1
3428 6.8 0.28 0.43 7.6 0.030 30.0 110.0 0.99164 3.08 0.59 12.5 8 1
699 6.9 0.23 0.34 2.7 0.032 24.0 121.0 0.99020 3.14 0.38 12.4 7 1
830 7.4 0.32 0.22 1.7 0.051 50.0 179.0 0.99550 3.28 0.69 8.9 5 1
1722 7.1 0.16 0.44 2.5 0.068 17.0 31.0 0.99328 3.35 0.54 12.4 6 0
1948 8.8 0.35 0.49 1.0 0.036 14.0 56.0 0.99200 2.96 0.33 10.5 4 1
3230 6.7 0.19 0.41 15.6 0.056 75.0 155.0 0.99950 3.20 0.44 8.8 6 1
3563 7.2 0.55 0.09 1.5 0.108 16.0 151.0 0.99380 3.07 0.57 9.2 4 1
In [31]:
test[:10]
Out[31]:
features label1 label2
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
4530 6.4 0.240 0.32 14.90 0.047 54.0 162.0 0.99680 3.28 0.50 10.2 6 1
1627 5.6 0.120 0.26 4.30 0.038 18.0 97.0 0.99477 3.36 0.46 9.2 5 1
801 7.0 0.400 0.25 1.80 0.050 51.0 189.0 0.99174 3.00 0.55 11.4 6 1
4045 6.8 0.120 0.31 5.20 0.045 29.0 120.0 0.99420 3.41 0.46 9.8 7 1
3651 7.9 0.460 0.40 10.10 0.168 19.0 184.0 0.99782 3.06 0.62 9.5 5 1
2649 6.8 0.260 0.43 11.75 0.045 53.0 198.0 0.99690 3.26 0.55 9.5 5 1
3930 7.5 0.570 0.02 2.60 0.077 11.0 35.0 0.99557 3.36 0.62 10.8 6 0
3287 11.9 0.380 0.51 2.00 0.121 7.0 20.0 0.99960 3.24 0.76 10.4 6 0
5123 7.1 0.250 0.28 1.20 0.040 31.0 111.0 0.99174 3.18 0.53 11.1 5 1
3195 8.3 0.615 0.22 2.60 0.087 6.0 19.0 0.99820 3.26 0.61 9.3 5 0
In [32]:
print("Training set (filas,columnas): " + str(train.shape))
print("Test set (filas,columnas): " + str(test.shape))
Training set (filas,columnas): (3990, 13)
Test set (filas,columnas): (1330, 13)

Ahora que el conjunto de test ya está generado, no lo tenemos que tocar más hasta el final, una vez dispongamos de un modelo. Únicamente trabajaremos con el conjunto de train, por lo que es momento de decidir si hay que modificar nuestro conjunto de train, es decir, si habría que eliminar outliers, descartar features por elevada correlación, estandarizarlas o bien, seleccionar features a través de la reducción de dimensionalidad.

Las conclusiones al respecto son las siguientes:

  • No vale la pena seleccionar features o realizar reducción de dimensionalidad ya que disponemos solamente de 11 variables dependientes más las 2 independientes o objetivo, de todas formas, lo añadiremos en nuestros modelos para que seleccionen automáticamente las mejores variables y así, observaremos las variaciones en los resultados.

    La selección de features se suele realizar cuando son datasets con muchísimas features, utilizando normalmente los árboles decisión, ya que estos automáticamente te indican las mejores variables para el modelo. Si un árbol decide no hacer split por una determinada feature nunca, es una señal de que dicha feature probablemente sea irrelevante para nuestro problema.

    Por otro lado la reducción de dimensionalidad lo que hace son pliegues de las variables en la dirección de la máxima varianza posible. Por lo que nos genera menos features, según las que nosotros solicitemos como resultantes (PCA - Principal Component Analysis). Genera features sintéticas a partir de unas originales, así se genera el top principal component del dataset original. Esto es muy utilizado cuando disponemos de datasets con muchísimas features, ya que no podemos tratar con todas en nuestro modelo, distorsionarian la predicción del mismo.

  • En cuanto a la estandarización, se debe realizar siempre en los modelos de sklearn a excepción de los árboles (y parecidos, como Random Forest y GBTs) y en Naïve Bayes. Por lo que en el momento de probar los modelos, deberemos añadir el StandarScaler o no.
  • Como hemos visto con los histogramas, no existen outliers significativos que puedan distorsionar nuestro modelo. No eliminaremos ningún dato del dataset.
  • Por último, como no disponemos de correlaciones entre pares de variables muy elevadas, descartamos el eliminar variables por alta correlacion.
In [33]:
# vino_data será un dataframe obtenido mediante el método copy()
vino_data = train.copy() 

Para acabar de afirmar que no existen outliers significativos o relevantes en grandes cantidades, analizamos la variable "density", que es una de las más relevantes de nuestro dataset dada su correlación con la resta. Observamos que apenas hay una observación alejada, por lo que no es necesario eliminarla, ya que no conocemos si la misma puede afectar negativa o positivamente a nuestro modelo.

In [34]:
from scipy.cluster.vq import kmeans, vq

centroids, avg_distance = kmeans(train["features", "density"], 5)
groups, cdist = vq(train["features","density"], centroids)
In [35]:
y = np.arange(0,train["features","density"].shape[0])
plt.scatter(train["features","density"],  y , c=groups)
plt.xlabel("Grado")
plt.ylabel("Observaciones")
plt.show()

Modelos clasificadores

Los problemas clasificadores se tratan para solucionar cuestiones en las que necesitamos "etiquetar" dentro de una serie de grupos controlados los datos, para hacer las cosas aún más sencillas.

El objetivo es predecir a través de unas variables independientes o features, la variable label de cada observación. Por lo que si tenemos todos los features, podremos predecir el label de cada observación.

Creación de los Pipelines

Los Pipelines son la manera más sencilla de encadenar operaciones, es decir, a través de un Pipeline podemos configurar que por ejemplo haga la selección de features, una vez hecha, haga un árbol de decisión para nuestro modelo. Todo esto se encadena gracias a un Pipeline.

La finalidad de los mismos es agilizar el proceso y no tener que hacerlo todo paso a paso a mano, con el objetivo de simplificar.

Podemos generar multitud de combinaciones en nuestros modelos a la hora de realizar los Pipelines, como puede ser añadir o no la selección de features, que la haga con una función u otra, etc.

In [36]:
from sklearn.pipeline import Pipeline

1) Partiremos con los árboles de decisión:

In [37]:
# No es necesario meter en una Pipeline ya que realiza 
# automáticamente los pasos por si mismo.
from sklearn.tree import DecisionTreeClassifier

arbol = DecisionTreeClassifier()

2) Seguido del RandomForest:

Provinien del Ensemble Learning, es decir, consiste en combinar varios modelos en uno solo, con el fin de que el modelo final sea más robusto y generalice mejor.

Mediante varios modelos a la vez, los fallos de uno de los modelos son compensados por los aciertos de otros de los modelos sobre los mismos puntos.

Hay distintias técnicas de ensembles:

  • Bagging: en este caso se lanzan varios modelos a la vez (normalmente cada uno de ellos con un subconjunto del training set). Las predicciones finales se basan en el consenso al que llegan los modelos.
  • Boosting: en este caso se lanza primero un modelo. Posteriormente, se lanza otro utilizando la parte de las observaciones bien clasificadas por el primero, y las observaciones mal clasificadas por éste. Así sucesivamente. Por entenderlo mejor, el segundo modelo intenta corregir los fallos del primero, así sucesivamente.

En nuestro caso probaremos tanto con la técnica de Random Forest (Bagging de árboles de decisión) como con la Gradient BoostedTrees (Boosting de árboles), por tener algo más de variedad.

En este caso no estandarizaremos ni haremos selección de features, lo hará el mismo modelo automáticamente.

In [38]:
from sklearn.ensemble import RandomForestClassifier

random_forest = RandomForestClassifier()
In [39]:
from sklearn.ensemble import GradientBoostingClassifier

gradient_boosting = GradientBoostingClassifier()

3) Ahora vamos con las regresiones logísticas:

Es un algorítmo de clasificación, ya que clasifica en números entre 0 y 1. Gracias a esto, podemos interpretar el resultado de la predicción como un resultado probabilístico. Por lo general si el resultado >0.5 se predice como 1, y si el resultado es <0.5 se predice como 0.

En este caso destacaremos la estandarización y la selección de features:

In [40]:
# Iniciamos con una regresión logística, haremos dos pipelines, cada uno dispondrá de 
# un seleccionador de variables distinto, aparte de la correspondiente estandarización.

from sklearn.preprocessing import StandardScaler
from sklearn.feature_selection import RFECV, SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression

logreg_rfecv = Pipeline(steps=[("scaler",StandardScaler()),
                               # método de selección de features
                               ("rfecv",RFECV(estimator=LogisticRegression())), 
                               ("logreg",LogisticRegression())
                              ]
                       )

logreg_kbest = Pipeline(steps=[("scaler",StandardScaler()),
                               ("kbest",SelectKBest()),  # método de selección de features
                               ("logreg",LogisticRegression())
                              ]
                       )

4) Seguimos con Nearest Neighbors:

Es un modelo no lineal, el concepto en el que se basa es el siguiente: si conocemos el label/target/objetivo a predecir de un determinado punto, los puntos cercanos deberían tener el mismo label; esa será su predicción.

El resultado dependerá del número de vecinos que contemplemos para un punto a clasificar. Por lo que si en nuestro ejemplo, hay que considerar los 10 puntos más cercanos al que queremos clasificar y 8 de ellos son vino blanco, se predicirá que el punto a clasificar es vino blanco.

In [41]:
# Utilizamos Nearest Neighbors: con selección KBest y, sin selección:

from sklearn.neighbors import KNeighborsClassifier

neighbors_kbest = Pipeline(steps=[("scaler",StandardScaler()),
                                  ("kbest",SelectKBest()),
                                  ("knn",KNeighborsClassifier())
                                ]
                         )

neighbors = Pipeline(steps=[("scaler",StandardScaler()),
                            ("knn",KNeighborsClassifier())
                          ]
                   )

5) A continuación los SVMs o Support Vector Machines:

Son globalmente conocidos como de los mejores clasificadores que hay, por su elevado grado y nivel matemático.

In [42]:
# Utilizamos SVMs sin selección y con KBest:

from sklearn.svm import SVC

svm = Pipeline(steps=[("scaler",StandardScaler()),
                      ("svm",SVC())
                     ]
              )

svm_kbest = Pipeline(steps=[("scaler",StandardScaler()),
                            ("kbest",SelectKBest()),
                            ("svm",SVC())
                           ]
                    )

6) En penúltimo lugar, vamos con Naïve Bayes:

Son modelos puramente probabilísticos, muy utilizados en problemas de clasificación. Su teorema se basa en: pensar que todas las features son independientes entre si. De ahí que sea puramente probabilistico y no comprenda la verosimilitud como los otros modelos.

In [43]:
# Utilizamos Naïve Bayes con y sin KBest (en este caso no estandarizamos): 

from sklearn.naive_bayes import GaussianNB

nb = GaussianNB()

nb_kbest = Pipeline(steps=[("kbest",SelectKBest()),
                           ("nb",GaussianNB())])

7) Concluímos con Percepción Multicapa:

Es la red neuronal más sencilla que xiste, la manera de entenderlo a modo simple es como una sucesión de capas de regresiones logísticas.

In [44]:
# La percepción multicapa, sin selección:

from sklearn.neural_network import MLPClassifier

mlp = Pipeline(steps=[("scaler",StandardScaler()),
                      ("mlp",MLPClassifier())
                     ]
              )

Una vez tenemos todos los Pipelines descritos, habrá que seleccionar hiperparámetros a probar para cada uno de ellos, es decir, seleccionar hasta donde queremos que nuestros modelos entrenen y darles una guia. Por lo que habrá que indicar la profundidad de los árboles de decisió, en cuantas features queremos que se quede cierto modelo, etc.

Selección de hiperparámetros

Para el árbol de decisión: decidimos entre profundidades del 1 al 10.

In [45]:
grid_arbol = {"max_depth":list(range(1,11))
             }

Para el random forest: Se suelen no necesitar más de 150 árboles, el modelo no empeora por exceso de árboles, pero apenas gana rendimiento. En este caso probamos distintas profundidades al azar.

In [46]:
grid_random_forest = {"n_estimators": [150],
                      "max_depth": [3,6,12,17,24,30],
                      "max_features": ["sqrt", 3,4]
                     }

Para el Gradient Boosting: se suele utilizar por norma general "deviance" en la función "loss". En cuanto al learning_rate la norma es; cuanto más alta más aportará cada nuevo árbol, siempre dependiendo del número de árboles que añadamos, por lo que un learning_rate alto junto a un n_estimators alta, puede dar sobreajuste. Si ponemos un número elevado de árboles podría sobreajustar, no como en Random Forest.

En boosting, a no ser que dispongamos de muchas features, los árboles deben tener poca profundidad, esto se hace así con el objetivo de que cada pequeño árbol generado se vaya rectificando con otro pequeño árbol de manera secuencial, por lo que la colección final de pequeños árboles consigue captar muy bien la tendencia de los datos.

Para "subsample" se conoce lo siguiente: para valores por debajo de 1.0 (<1.0), cada nuevo árbol entrenará solo con el porcentaje del training set especificado; esto puede ayudar a evitar la sobre-especialización de los árboles y por tanto, bajar nuestro variance del modelo.

Por concluir, el "max_features" funciona de igual manera que en el Random Forest.

In [47]:
grid_gradient_boosting = {"loss": ["deviance"],
                          "learning_rate": [0.05, 0.2, 0.6],
                          "n_estimators": [15, 35, 110, 180],
                          "max_depth": [1, 2, 3, 4, 5],
                          "subsample": [1.0, 0.7, 0.4],
                          "max_features": ["sqrt", 3, 4]
                         }

De igual manera, procedemos con la selección de hiperparámetros de las regresiones logísticas:

El primer paso (rfecv_step) lo utilizamos para indicar que queremos que pruebe a quitar features de una en una, es la manera más conservadora de realizarlo. En cambio, con "rfecv__cv" estamos indicando que queremos que realize 6 pliegues (folds) en la CV interna del RFECV.

Indicamos las regularizaciones y su fuerza con "logregpenalty" y "logregC".

Tambien debemos indicar con "logreg__max_iter" el número de iteraciones que queremos que realice, no suele preocupar si es un número elevado en el caso de regresiones lineales y logísticas.

Para ejecutarlo seleccionamos "liblinear" en "logreg__solver" ya que suele ser el método más rápido.

In [48]:
grid_logreg_rfecv = {"rfecv__step": [1],
                     "rfecv__cv": [6],
                     "logreg__penalty": ["l1", "l2"],
                     "logreg__C": [0.1, 0.9, 2.0, 6.0],
                     "logreg__fit_intercept": [True],
                     "logreg__max_iter": [40, 130, 400],
                     "logreg__solver": ["liblinear"]
                    }

En este caso debemos indicar la manera de clasificar con kbest las features y, escoger el número "kbest__k" de features con las que quedarse, en nuestro modelo indicaremos que se quede con 7 (no conviene sobre eliminar features des de mi punto de vista cuando tampoco tenemos tantas).

La resta de pasos funcionan como en la regresión anterior, explicados con el modelo de regresión logística y selección de features RFECV.

In [49]:
grid_logreg_kbest = {"kbest__score_func": [f_classif],
                     "kbest__k": [1,2,3,4,5,6,7],
                     "logreg__penalty": ["l1", "l2"],
                     "logreg__C": [0.1, 0.9, 2.0, 6.0],
                     "logreg__fit_intercept": [True],
                     "logreg__max_iter": [40, 130, 400],
                     "logreg__solver": ["liblinear"]
                    }

Seguimos con los Nearest Neighbors: otro modelo de clasificación categórica, al cual hay que pasarle valores "k" positivos impares para que conozco la mejor predicción. No es lineal

En el caso del modelo con selección kbest indicamos que queremos quedarnos con 6 features.

Por otra banda, solo deberemos seleccionar el numero de "k" que queremos incorporarle, y indicarle con "weights" la manera de ponderar o no las clasificaciones en función de la inversa de la distancia a cada vecino.

In [50]:
grid_neighbors_kbest = {"kbest__score_func": [f_classif],
                        "kbest__k": [1,2,3,4,5,6],
                        "knn__n_neighbors": [3, 5, 7, 9, 11, 13, 15],
                        "knn__weights": ["uniform","distance"]
                       }
In [51]:
grid_neighbors = {"knn__n_neighbors": [3, 5, 7, 9, 11, 13, 15],
                  "knn__weights": ["uniform","distance"]
                 }

Ahora realizaremos el modelo Support Vector Machines (SVM): el cual lo que procura es realizar en el caso de que el clasificador sea lineal, realizar una división tal que se maximice el margen entre los dos tipos de clases, con esto intenta dejar equitativamente la división entre categorías.

En el caso de ser un modelo no lineal, cogerá el dataset y le aplicará una transformación no lineal, llevandolo a una dimensión espacial superior. En una de estas dimensiones se podrá realizar una división lineal, que una vez encontrada, el modelo transformará nuevamente a una no lineal ya con la división (clasificación por categorías) realizada de manera automática.

Representamos los pasos que sigue el modelo visualmente:

Seleccionaremos los hiperparámetros para el modelo SVM sin selección de features y con seleccion de features siguiendo KBest.

In [52]:
grid_svm = {"svm__C": [0.01, 0.1, 0.5, 1.0, 5.0, 20.0],
            "svm__kernel": ["linear"],      # No calculamos polinómicos ni gaussianos: "poly","rbf"
            "svm__degree": [2,3,4,5],       # Por el tiempo de ejecución
            "svm__gamma": [0.001, 0.1, "auto", 1.0, 10.0, 20.0]
           }

grid_svm_kbest = {"kbest__score_func": [f_classif],
                  "kbest__k": [1,2,3,4,5,6,7,8,9],        # Indicamos que queden 9 features
                  "svm__C": [0.01, 0.1, 0.5, 1.0, 5.0, 20.0],
                  "svm__kernel": ["linear"],  # No calculamos polinómicos ni gaussianos: "poly","rbf".                
                  "svm__degree": [2,3,4,5],   # Por el tiempo de ejecución
                  "svm__gamma": [0.001, 0.1, "auto", 1.0, 10.0, 20.0]
                 }

Ya casi estamos en cuanto a la seleccion de hiperparámetros para nuestros modelos, ahora vamos con Naïve Bayes: su mayor ventaja es que entrena y predice muy rápido en comparación a la resta, su tiempo de ejecución es corto. En vez de tratar la probabilidad de un registro, va más allá, trata de la probabilidad de que haya la probabilidad de un registro. No es lineal, sino cuadrático.

Como hemos realizado los pipelines con y sin selección de KBest, solo sera necesario configurar los hiperparámetros de KBest, por lo que nos decantaremos por dejar 9 features.

In [53]:
grid_nb_kbest = {"kbest__score_func": [f_classif],
                 "kbest__k": [1,2,3,4,5,6,7,8,9]
                }

Para concluir, acabaremos con el modelo de percepción multicapa o MLP: se entiende como deep learning o redes neuronales. El objetivo del modelo es crear pequeñas neuronas que, van realizando predicciones mejores que las precedesoras hasta obtener un resultado.

Vamos a introducir los hiperparámetros:

A través de "mlp__hidden_layer_sizes" seleccionamos distintas arquitecturas de la red en las que puede ir trabajando hasta obtener un resultado, podríamos poner infinitas si lo deseamos.

Decidimos activar con "mlp__activation" la función logística, la función de activación lineal rectificada (ReLu) para redes de alto rendimiento y, la función tanh, está última contiene una serie de funciones que se utilizan para operaciones matemáticas.

En el caso de las redes neuronales, el mejor método para "solver" suele ser "adam".

No pondremos mucha regularización en las redes con "alpha". En cuanto al "mlp__validation_fraction"; este se utiliza para indicar el número de iteraciones de Gradient Descent, muchas de estas pueden llegar a sobreajuste, por lo que siempre es mejor no esperar a que converja el gradient descent (conocido como Early Stopping).

Junto con el "validation_fraction", indicaremos "True" en el "early_stopping" tal y como hemos comentado. Como hemos indicado que pare antes de que converja, no importa el "max_iter" que indiquemos ya que probablemente pare antes, así que pondremos uno elevado.

Por último le indicamos con "learning_rate_init" el tamaño de los pasos y el learning rate inicial.

In [54]:
grid_mlp = {"mlp__hidden_layer_sizes": [(6,),
                                        (6,6),
                                        (70,),
                                        (70,70),
                                        (130,),
                                        (130,130,130)],
            "mlp__activation": ["logistic", "relu", "tanh"],
            "mlp__solver": ["adam"],
            "mlp__alpha": [0.0, 0.001, 0.1],
            "mlp__validation_fraction": [0.1],
            "mlp__early_stopping": [True],
            "mlp__max_iter": [5500],
            "mlp__learning_rate_init": [0.001, 0.2, 0.7]
           }

Grid Search y selección del modelo

De todos los hiperparámetros propuestos a nuestros modelos anteriormente, el Grid Search es la técnica más sencilla para buscar las mejores combinaciones de los mismos para cierto modelo y dataset.

Principalmente, genera automáticamente todas las combinaciones de todos los valores que le hemos dicho que pruebe (es lo conocido como grid), y así, mira qué resultados obtiene cada una de las validaciones (cross validation).

Realizamos el código correspondiente al Grid Search para cada modelo diseñado:

Esto lo haremos construyendo cada GridSearchCV con su Pipeline (configuradas previamente) y su grid de hiperparámetros (configurados previamente).

Lo que haremos sera indicarle a Grid Search que elija el modelo con mayor valor de accuracy.

In [55]:
from sklearn.model_selection import GridSearchCV

gs_arbol = GridSearchCV(arbol,
                       grid_arbol,
                       cv=10,
                       scoring="accuracy",
                       verbose=1,
                       n_jobs=-1)

gs_random_forest = GridSearchCV(random_forest,
                                grid_random_forest,
                                cv=10,
                                scoring="accuracy",
                                verbose=1,
                                n_jobs=-1)

gs_gradient_boosting = GridSearchCV(gradient_boosting,
                                    grid_gradient_boosting,
                                    cv=10,
                                    scoring="accuracy",
                                    verbose=1,
                                    n_jobs=-1)

gs_logreg_rfecv = GridSearchCV(logreg_rfecv,
                               grid_logreg_rfecv,
                               cv=10,
                               scoring="accuracy",
                               verbose=1,
                               n_jobs=-1)

gs_logreg_kbest = GridSearchCV(logreg_kbest,
                               grid_logreg_kbest,
                               cv=10,
                               scoring="accuracy",
                               verbose=1,
                               n_jobs=-1)

gs_neighbors = GridSearchCV(neighbors,
                            grid_neighbors,
                            cv=10,
                            scoring="accuracy",
                            verbose=1,
                            n_jobs=-1)

gs_neighbors_kbest = GridSearchCV(neighbors_kbest,
                                  grid_neighbors_kbest,
                                  cv=10,
                                  scoring="accuracy",
                                  verbose=1,
                                  n_jobs=-1)

gs_svm = GridSearchCV(svm,
                      grid_svm,
                      cv=10,
                      scoring="accuracy",
                      verbose=1,
                      n_jobs=-1)

gs_svm_kbest = GridSearchCV(svm_kbest,
                            grid_svm_kbest,
                            cv=10,
                            scoring="accuracy",
                            verbose=1,
                            n_jobs=-1)

gs_nb = GridSearchCV(nb,
                     {},  # No hay grid
                     cv=10,
                     scoring="accuracy",
                     verbose=1,
                     n_jobs=-1)
               
gs_nb_kbest = GridSearchCV(nb_kbest,
                           grid_nb_kbest,
                           cv=10,
                           scoring="accuracy",
                           verbose=1,
                           n_jobs=-1)

gs_mlp = GridSearchCV(mlp,
                      grid_mlp,
                      cv=10,
                      scoring="accuracy",
                      verbose=1,
                      n_jobs=-1)

Una vez está todo realizado, procedemos a meter todos los GridSearchCV dentro de un diccionario, para así poder ver los resultados de una manera mucho más práctica e intuitiva.

In [56]:
all_grid_searchs = {"gs_arbol":gs_arbol,
                    "gs_random_forest":gs_random_forest,
                    "gs_gradient_boosting":gs_gradient_boosting,
                    "gs_logreg_rfecv":gs_logreg_rfecv,
                    "gs_logreg_kbest":gs_logreg_kbest,
                    "gs_neighbors":gs_neighbors,
                    "gs_neighbors_kbest":gs_neighbors_kbest,
                    "gs_svm":gs_svm,
                    "gs_svm_kbest":gs_svm_kbest,
                    "gs_nb":gs_nb,
                    "gs_nb_kbest":gs_nb_kbest,
                    "gs_mlp":gs_mlp}

Una vez guardados en un diccionario, por lo que a través de la iteración con un bucle podremos ejecutar todos los Grid Search de cada uno de los modelos generados. Y para cada uno de ellos, encontraremos la mejor combinación de hiperparámetros.

In [57]:
for nombre, grid_search in all_grid_searchs.items():
    print("Haciendo Grid Search de %s..." % nombre)
    grid_search.fit(train["features"], train["label2", "color"])
Haciendo Grid Search de gs_arbol...
Fitting 10 folds for each of 10 candidates, totalling 100 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 100 out of 100 | elapsed:    0.9s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
Haciendo Grid Search de gs_random_forest...
Fitting 10 folds for each of 18 candidates, totalling 180 fits
[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed:    8.1s
[Parallel(n_jobs=-1)]: Done 180 out of 180 | elapsed:   46.9s finished
Haciendo Grid Search de gs_gradient_boosting...
Fitting 10 folds for each of 540 candidates, totalling 5400 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 125 tasks      | elapsed:    5.6s
[Parallel(n_jobs=-1)]: Done 412 tasks      | elapsed:   19.1s
[Parallel(n_jobs=-1)]: Done 852 tasks      | elapsed:   51.6s
[Parallel(n_jobs=-1)]: Done 1356 tasks      | elapsed:  1.6min
[Parallel(n_jobs=-1)]: Done 1840 tasks      | elapsed:  2.7min
[Parallel(n_jobs=-1)]: Done 2762 tasks      | elapsed:  3.7min
[Parallel(n_jobs=-1)]: Done 3582 tasks      | elapsed:  5.2min
[Parallel(n_jobs=-1)]: Done 4808 tasks      | elapsed:  6.7min
[Parallel(n_jobs=-1)]: Done 5400 out of 5400 | elapsed:  8.0min finished
Haciendo Grid Search de gs_logreg_rfecv...
Fitting 10 folds for each of 24 candidates, totalling 240 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed:   11.9s
[Parallel(n_jobs=-1)]: Done 192 tasks      | elapsed:   52.9s
[Parallel(n_jobs=-1)]: Done 240 out of 240 | elapsed:  1.1min finished
Haciendo Grid Search de gs_logreg_kbest...
Fitting 10 folds for each of 168 candidates, totalling 1680 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 200 tasks      | elapsed:    1.3s
[Parallel(n_jobs=-1)]: Done 1400 tasks      | elapsed:    8.7s
[Parallel(n_jobs=-1)]: Done 1680 out of 1680 | elapsed:   10.7s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
Haciendo Grid Search de gs_neighbors...
Fitting 10 folds for each of 14 candidates, totalling 140 fits
[Parallel(n_jobs=-1)]: Done  60 tasks      | elapsed:    1.3s
[Parallel(n_jobs=-1)]: Done 140 out of 140 | elapsed:    2.8s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
Haciendo Grid Search de gs_neighbors_kbest...
Fitting 10 folds for each of 84 candidates, totalling 840 fits
[Parallel(n_jobs=-1)]: Done 128 tasks      | elapsed:    1.1s
[Parallel(n_jobs=-1)]: Done 728 tasks      | elapsed:    7.5s
[Parallel(n_jobs=-1)]: Done 840 out of 840 | elapsed:    8.9s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
Haciendo Grid Search de gs_svm...
Fitting 10 folds for each of 144 candidates, totalling 1440 fits
[Parallel(n_jobs=-1)]: Done  76 tasks      | elapsed:    1.3s
[Parallel(n_jobs=-1)]: Done 376 tasks      | elapsed:    6.0s
[Parallel(n_jobs=-1)]: Done 876 tasks      | elapsed:   14.1s
[Parallel(n_jobs=-1)]: Done 1440 out of 1440 | elapsed:   38.1s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
Haciendo Grid Search de gs_svm_kbest...
Fitting 10 folds for each of 1296 candidates, totalling 12960 fits
[Parallel(n_jobs=-1)]: Done  76 tasks      | elapsed:    3.0s
[Parallel(n_jobs=-1)]: Done 376 tasks      | elapsed:   13.2s
[Parallel(n_jobs=-1)]: Done 876 tasks      | elapsed:   26.4s
[Parallel(n_jobs=-1)]: Done 1576 tasks      | elapsed:   47.4s
[Parallel(n_jobs=-1)]: Done 2476 tasks      | elapsed:  1.1min
[Parallel(n_jobs=-1)]: Done 3576 tasks      | elapsed:  1.5min
[Parallel(n_jobs=-1)]: Done 4876 tasks      | elapsed:  2.1min
[Parallel(n_jobs=-1)]: Done 6376 tasks      | elapsed:  2.7min
[Parallel(n_jobs=-1)]: Done 8076 tasks      | elapsed:  3.4min
[Parallel(n_jobs=-1)]: Done 9976 tasks      | elapsed:  4.3min
[Parallel(n_jobs=-1)]: Done 12168 tasks      | elapsed:  5.1min
[Parallel(n_jobs=-1)]: Done 12960 out of 12960 | elapsed:  5.4min finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
Haciendo Grid Search de gs_nb...
Fitting 10 folds for each of 1 candidates, totalling 10 fits
Haciendo Grid Search de gs_nb_kbest...
Fitting 10 folds for each of 9 candidates, totalling 90 fits
[Parallel(n_jobs=-1)]: Done  10 out of  10 | elapsed:    0.0s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done  90 out of  90 | elapsed:    0.4s finished
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
Haciendo Grid Search de gs_mlp...
Fitting 10 folds for each of 162 candidates, totalling 1620 fits
[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed:    3.6s
[Parallel(n_jobs=-1)]: Done 192 tasks      | elapsed:   40.6s
[Parallel(n_jobs=-1)]: Done 442 tasks      | elapsed:  1.6min
[Parallel(n_jobs=-1)]: Done 792 tasks      | elapsed:  3.2min
[Parallel(n_jobs=-1)]: Done 1242 tasks      | elapsed:  5.1min
[Parallel(n_jobs=-1)]: Done 1620 out of 1620 | elapsed:  7.0min finished

Una vez finalizado, vamos a ver cuál ha sido el mejor resultado de la cross validation en cada GridSearchCV realizado, lo haremos de la siguiente manera:

In [59]:
best_score_gridsearch = [(modelo, grid_search.best_score_)
                         for modelo, grid_search
                         in all_grid_searchs.items()]

best_score_gridsearch
Out[59]:
[('gs_arbol', 0.9824561403508772),
 ('gs_random_forest', 0.9932330827067668),
 ('gs_gradient_boosting', 0.9942355889724309),
 ('gs_logreg_rfecv', 0.9927318295739347),
 ('gs_logreg_kbest', 0.987719298245614),
 ('gs_neighbors', 0.9924812030075186),
 ('gs_neighbors_kbest', 0.9892230576441101),
 ('gs_svm', 0.9937343358395989),
 ('gs_svm_kbest', 0.9919799498746867),
 ('gs_nb', 0.968671679197995),
 ('gs_nb_kbest', 0.9691729323308271),
 ('gs_mlp', 0.9949874686716791)]

Vamos a introduciros en un DataFrame de pandas para que sea más senzillo de visualizar y trastear:

In [60]:
best_score_gridsearch_df = pd.DataFrame(best_score_gridsearch,
                                        columns=["GridSearchCV", "Mejor Score"])
# Lo vamos a poner en orden

best_score_gridsearch_df_sort = (best_score_gridsearch_df.sort_values(by="Mejor Score", ascending=False))

best_score_gridsearch_df_sort
Out[60]:
GridSearchCV Mejor Score
11 gs_mlp 0.994987
2 gs_gradient_boosting 0.994236
7 gs_svm 0.993734
1 gs_random_forest 0.993233
3 gs_logreg_rfecv 0.992732
5 gs_neighbors 0.992481
8 gs_svm_kbest 0.991980
6 gs_neighbors_kbest 0.989223
4 gs_logreg_kbest 0.987719
0 gs_arbol 0.982456
10 gs_nb_kbest 0.969173
9 gs_nb 0.968672
In [61]:
print("El mejor modelo seleccionado por GridSearchCV es: " + 
      str(best_score_gridsearch_df_sort["GridSearchCV"].iloc[0]) + 
      " con un score de " + str((best_score_gridsearch_df_sort["Mejor Score"].iloc[0])*100) + 
      "%")
El mejor modelo seleccionado por GridSearchCV es: gs_mlp con un score de 99.49874686716791%

La percepción multicapa es el modelo vencedor, el cual es el mejor grid search en validación cruzada. El score nos devuelve un número entre 0 y 1 (lo hemos múltiplicado por 100 para tener el valor en porcentaje), el cual nos indica el porcentaje de observacions bien clasificadas. Por otra banda, todos los modelos disponen de scores muy elevadas, pudiendo afirmar que hay 7 modelos por encima del 99% de score, unas cifras muy buenas en este caso para el problema categórico.

Por lo que escogemos el mejor GridSearchCV, que en este caso es el GradientBoosting:

In [62]:
best_gridsearch = all_grid_searchs["gs_mlp"]

Ahora dentro de este GridSearchCV, el mejor PipeLine (modelo) ha sido:

In [63]:
best_pipeline = best_gridsearch.best_estimator_

best_pipeline
Out[63]:
Pipeline(steps=[('scaler', StandardScaler()),
                ('mlp',
                 MLPClassifier(activation='logistic', alpha=0.001,
                               early_stopping=True, hidden_layer_sizes=(6, 6),
                               learning_rate_init=0.2, max_iter=5500))])

Una vez disponemos del modelo seleccionado, aúnque Scikit-Learn entrena todo el conjunto de train con el modelo que mejor valoración tiene del GridSearchCV, debido al argumento opcional refit=True que viene marcado por defecto, vamos a proceder a entrenarlo nuevamente para no dejar duda al respecto.

In [64]:
best_pipeline.fit(train["features"], train["label2","color"])
Out[64]:
Pipeline(steps=[('scaler', StandardScaler()),
                ('mlp',
                 MLPClassifier(activation='logistic', alpha=0.001,
                               early_stopping=True, hidden_layer_sizes=(6, 6),
                               learning_rate_init=0.2, max_iter=5500))])

Mediciones de la predicción sobre el conjunto de test

Primero realizaremos las métricas, de las que obtenemos resultados númericos:

Accuracy:

Vamos a extraer el porcentaje de observacions que el modelo clasifica correctamente en el conjunto de test:

In [71]:
from sklearn.metrics import accuracy_score

accuracy_en_test = accuracy_score(y_true = test["label2", "color"],
                                  y_pred = best_pipeline.predict(test["features"])
                                 )

print("Nuestro modelo tiene un accuracy en el conjunto de test de %s " % accuracy_en_test)
Nuestro modelo tiene un accuracy en el conjunto de test de 0.9932330827067669 

De todas formas, no siempre el accuracy/score es la mejor manera de clasificar, ya que no discrimina y diferencia entre lo que debería ser un uno de lo que debería ser un cero.

Existen otras métricas para valorar el rendimiento de nuestro clasificador:

La matriz de confusión:

Se utiliza normalmente en problemas de clasificación binaria como nuestro caso. Se trata de una tabla con cuatro cuadrantes en nuestro caso:

- El cuadrante superior izquierdo muestra el número de **unos** bien clasificados.
- El cuadrante inferior izquierdo muestra el número de *ceros* clasificados como *uno*.
- El cuadrante superior derecho muestra el número de *unos* clasificados como *cero*.
- El cuadrante inferor derecho muestra el número de *ceros* bien clasificados.

Vamos a proceder a codificar la matriz de confusión gracias a la función de sklearn.metrics:

In [72]:
from sklearn.metrics import confusion_matrix

matriz_confusion = confusion_matrix(y_true = test["label2", "color"],
                                    y_pred = best_pipeline.predict(test["features"])
                                   )

matriz_confusion
Out[72]:
array([[337,   7],
       [  2, 984]], dtype=int64)

Podemos hacer la matriz de confusión algo más bonita, por lo que la pasaremos a un DataFrame de Pandas.

In [73]:
matriz_conf_df = pd.DataFrame(matriz_confusion)

# Cambiamos el nombre de los 0 y 1 a su valor original en formato string, Vino tinto y Vino blanco. (red, white)

matriz_conf_df = matriz_conf_df.rename(columns = {0: "Vino tinto", 1: "Vino blanco"})

matriz_conf_df.index = ["Vino tinto", "Vino blanco"]

# Indicamos lo que son columnas y filas:

matriz_conf_df.columns.name = "Predicho"
matriz_conf_df.index.name ="Real"

matriz_conf_df
Out[73]:
Predicho Vino tinto Vino blanco
Real
Vino tinto 337 7
Vino blanco 2 984

De manera gráfica, y ya pasamos a describir los resultados obtenidos:

In [74]:
plt.figure(figsize=(6,4))

sns.heatmap(matriz_conf_df,
            annot=True, 
            cmap="BuGn",
            fmt="d")
pass

Vamos a ver el valor relativo de cada tipo de vino sobre el total del conjunto test para poder evaluarlo posteriormente:

In [75]:
test["label2", "color"].value_counts()
Out[75]:
1    986
0    344
Name: (label2, color), dtype: int64

Podemos de este modo, dar una conclusión sobre nuestro modelo escogido (GradientBoosting) con GridSearchCV aplicado al conjunto de testing, siguiendo la matriz de confusión:

  • Ha predicho bien 337 vinos tintos (zeros), de un total de 344 -> Verdaderos Positivos.
  • Ha predicho 2 vino tinto como vino blanco -> Falso Positivo.
  • Ha predicho 7 vinos blancos como vinos tintos -> Falsos Negativos.
  • Ha predicho bien 984 vinos blancos (unos), de un total de 986 -> Verdaderos Negativos.

Apreciamos que la precisión en la predicción del modelo es muy elevada, su volumen de error es mínimo y no detectamos cantidades dispares en los errores predichos.

Precision:

Se entiende técnicamente como la habilidad que tiene un clasificador para no predecir los ceros como unos, es decir, el % de acierto de que los ceros sean ceros y los unos sean unos.

Vamos a ello:

In [76]:
from sklearn.metrics import precision_score

precision_best_pipeline = precision_score(y_true = test["label2", "color"],
                                    y_pred = best_pipeline.predict(test["features"])
                                   )

print("Precision: " + str(precision_best_pipeline))
print("Cuando nuestro clasificador predice una observación como uno, acierta el " + 
      str(precision_best_pipeline*100) + "% de los casos (en nuestro conjunto de testing.)")
Precision: 0.992936427850656
Cuando nuestro clasificador predice una observación como uno, acierta el 99.2936427850656% de los casos (en nuestro conjunto de testing.)

Recall:

Se entiende como la habilidad que tiene un clasificador para predecir todos los unos reales.

Lo ejecutamos de la siguiente manera:

In [77]:
from sklearn.metrics import recall_score

recall_best_pipeline = recall_score(y_true = test["label2", "color"],
                                    y_pred = best_pipeline.predict(test["features"])
                                   )

print("Recall: " + str(recall_best_pipeline))
print("Nuestro clasificador ha conseguido predecir bien el " + str(recall_best_pipeline) + 
      "% de los unos (en nuestro conjunto de testing).")
Recall: 0.9979716024340771
Nuestro clasificador ha conseguido predecir bien el 0.9979716024340771% de los unos (en nuestro conjunto de testing).

F1-Score:

Esta última métrica, lo que hace principalmente es combinar en una sola métrica la precision y recall a través de una fórmula matemática:

In [78]:
from sklearn.metrics import f1_score

f1_score_best_pipeline = f1_score(y_true = test["label2", "color"],
                                    y_pred = best_pipeline.predict(test["features"])
                                   )

print("F1-Score: " + str(f1_score_best_pipeline))
F1-Score: 0.9954476479514417

Dado el más que próximo valor a 1 de nuestro F1-Score, podemos decir que es una muy buena predicción, ya que significará que tambien el precision y el recall son muy cercanos a 1.

Por otra banda, tambien podemos determinar si un modelo es bueno a nivel de predicción via gráfica, por lo que vamos a ponernos manos a la obra con los distintos métodos gráficos que existen:

Curva ROC:

La curva ROC (Receiver Operating Characteristic) es un gráfico muy utilizado cuando queremos visualizar el rendimiento de un clasificador binario. Los ejes del gráfico son:

  • El eje x es el Falso Positivo Ratio (FPR)
  • El eje y es el Verdadero Positivo Ratio (TPR)

La manera matemática de interpretarlos queda tal que:

Los usos de la curva son principalmente estos:

  • Cuando los clasificadores generan scoring, determinamos que son ceros las observacions que puntúan por debajo de cierto valor umbral; y unos los que puntúan por encima. La curva ROC lo que nos permite es visualizar cuál es el mejor umbral que podemos poner. Por lo que según el umbral, unos quedaran por encima y otros por debajo, separando así la predicción.
  • Cuando los clasificadores no generan scoring, la curva ROC será solo una línea que une el (0,0), el punto (FPR, TPR) de nuestro modelo, y el (1,1). Por lo que medir el área bajo esta curva (área Under ROC) es una métrica de rendimiento que nos permite comparar clasificadores.

Nosotros vamos a centrarnos en el medir el área de abajo de la curva, por lo que utilizaremos el segundo uso.

Lo ideal es que un modelo tenga un False Positive Rate muy bajo (lo más cercano a 0 posible), y un True Positive Rate muy elevado (aproximado a 1). Si es así, el área bajo la curva ROC (AUC) será lo más elevada posible.

Vamos a dibujarla con nuestro modelo seleccionado (best_pipeline):

In [79]:
from sklearn.metrics import roc_curve # Hace automáticamente una tupla con 3 arrays de Numpy

roc = roc_curve(y_true = test["label2", "color"], 
                y_score = best_pipeline.predict(test["features"]),
                pos_label=1.0 # Indicamos que los casos positivos son los unos
                                   )
roc
Out[79]:
(array([0.        , 0.02034884, 1.        ]),
 array([0.       , 0.9979716, 1.       ]),
 array([2, 1, 0]))
In [80]:
# Cogemos las coordenadas del resultado, las extraemos del array generado:
coordenadas_fpr = roc[0]
coordenadas_tpr = roc[1]

plt.figure(figsize=(12,8))

# En primer lugar, pintamos el random guess:
plt.plot([0,1], [0,1], "g--")

# Ahora el FPR y TPR de nuestro modelo:
plt.plot(coordenadas_fpr, coordenadas_tpr, "r")

# Vamos a darle color al área bajo la curva:
plt.fill_between(coordenadas_fpr, coordenadas_tpr, alpha=0.3, color="orange")

plt.title("Curva ROC Mejor Modelo (Gradient Boosting)")
plt.xlabel("Falso Positivo Ratio")
plt.ylabel("Verdadero Positivo Ratio")
pass

Podemos observar como la curva ROC se despega muchísimo del random guess (linea diagonal discontinua) por arriba, siendo prácticamente una recta al inicia en su crecimiento.

Al tener la curva ROC muy por encima del rangom guess, ya estaríamos realizando una mejor predicción que si hiciesemos a boleo, y encima, en nuestro caso, rozamos el 1, por lo que nuestra predicción es prácticamente precisa al 100%.

Vamos a calcular el área (AUC), la cual debería ser prácticamente 1, ya que nuestra curva ROC ocupa todo el cuadro.

In [81]:
from sklearn.metrics import roc_auc_score

auc_best_pipeline = roc_auc_score(y_true = test["label2", "color"], 
                                  y_score = best_pipeline.predict(test["features"]))

print("AUC: " + str(auc_best_pipeline))
AUC: 0.9888113826123872

Finalmente así ha sido, el área bajo la curva ROC es prácticamente 1, rozando el máximo.

Modelos Regresores:

Los problemas de regresión, a diferencia de los de clasificación donde predecimos labels o clases, se utilizan para predecir un número; como siempre, en función de una serie de features que tiene cada observación obtenida en un dataset.

En nuestro caso, vamos a procurar predecir la calidad del vino (con valores del 1 al 10) de nuestro dataset a través de las 11 features de las que disponemos.

Por lo que seguiremos trabajando con el mismo dataset, pero en esta ocasión, no utilizaremos la variable a predecir en el caso de clasificación ("color"), simplemente nos quedaremos con las features que caracterizan y dan atributos a las observaciones relacionadas con distintos vinos y la actual label, quality.

In [82]:
df_vino1.head()
Out[82]:
features label1 label2
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality color
0 5.20 0.34 0.00 1.8 0.050 27.0 63.0 0.99160 3.68 0.79 14.0 6 0
1 6.20 0.55 0.45 12.0 0.049 27.0 186.0 0.99740 3.17 0.50 9.3 6 1
2 7.15 0.17 0.24 9.6 0.119 56.0 178.0 0.99578 3.15 0.44 10.2 6 1
3 6.70 0.64 0.23 2.1 0.080 11.0 119.0 0.99538 3.36 0.70 10.9 5 0
4 7.60 0.23 0.34 1.6 0.043 24.0 129.0 0.99305 3.12 0.70 10.4 5 1
In [83]:
df_vino_r = df_vino1.iloc[:, :-1]
In [84]:
df_vino_r.head()
Out[84]:
features label1
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality
0 5.20 0.34 0.00 1.8 0.050 27.0 63.0 0.99160 3.68 0.79 14.0 6
1 6.20 0.55 0.45 12.0 0.049 27.0 186.0 0.99740 3.17 0.50 9.3 6
2 7.15 0.17 0.24 9.6 0.119 56.0 178.0 0.99578 3.15 0.44 10.2 6
3 6.70 0.64 0.23 2.1 0.080 11.0 119.0 0.99538 3.36 0.70 10.9 5
4 7.60 0.23 0.34 1.6 0.043 24.0 129.0 0.99305 3.12 0.70 10.4 5
In [85]:
df_vino_r.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 5320 entries, 0 to 5319
Data columns (total 12 columns):
 #   Column                            Non-Null Count  Dtype  
---  ------                            --------------  -----  
 0   (features, fixed acidity)         5320 non-null   float64
 1   (features, volatile acitidy)      5320 non-null   float64
 2   (features, citric acid)           5320 non-null   float64
 3   (features, residual sugar)        5320 non-null   float64
 4   (features, chlorides)             5320 non-null   float64
 5   (features, free sulfur dioxide)   5320 non-null   float64
 6   (features, total sulfur dioxide)  5320 non-null   float64
 7   (features, density)               5320 non-null   float64
 8   (features, pH)                    5320 non-null   float64
 9   (features, sulphates)             5320 non-null   float64
 10  (features, alcohol)               5320 non-null   float64
 11  (label1, quality)                 5320 non-null   int64  
dtypes: float64(11), int64(1)
memory usage: 498.9 KB

Como ya habíamos limpiado el dataset anteriormente, no es necesario realizarlo nuevamente. Lo que haremos es volver a dividir nuestro dataset en train y test, para poder entrenar el conjunto de training y posteriormente aplicarlo a test para conocer nuestro porcentaje de precisión.

In [86]:
train_r, test_r = train_test_split(df_vino_r,
                                   train_size = 0.8,  # 80% training
                                   test_size = 0.2,   # 20% testing
                                   random_state=24) # seed para el generador
In [87]:
train_r[:10]
Out[87]:
features label1
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality
4118 6.4 0.31 0.40 6.4 0.039 39.0 191.0 0.99513 3.14 0.52 9.8 5
2923 7.7 0.43 0.28 4.5 0.046 33.0 102.0 0.99180 3.16 0.56 12.2 7
5013 11.6 0.42 0.53 3.3 0.105 33.0 98.0 1.00100 3.20 0.95 9.2 5
1385 6.1 0.46 0.32 6.2 0.053 10.0 94.0 0.99537 3.35 0.47 10.1 5
3265 8.4 0.31 0.29 3.1 0.194 14.0 26.0 0.99536 3.22 0.78 12.0 6
2498 6.9 0.38 0.38 13.1 0.112 14.0 94.0 0.99792 3.02 0.48 9.2 5
4599 5.1 0.23 0.18 1.0 0.053 13.0 99.0 0.98956 3.22 0.39 11.5 5
2698 6.8 0.26 0.56 11.9 0.043 64.0 226.0 0.99700 3.02 0.63 9.3 5
2823 6.1 0.55 0.15 9.8 0.031 19.0 125.0 0.99570 3.36 0.47 10.2 6
1419 7.2 0.24 0.29 2.2 0.037 37.0 102.0 0.99200 3.27 0.64 11.0 7
In [88]:
test_r[:10]
Out[88]:
features label1
fixed acidity volatile acitidy citric acid residual sugar chlorides free sulfur dioxide total sulfur dioxide density pH sulphates alcohol quality
4530 6.4 0.240 0.32 14.90 0.047 54.0 162.0 0.99680 3.28 0.50 10.2 6
1627 5.6 0.120 0.26 4.30 0.038 18.0 97.0 0.99477 3.36 0.46 9.2 5
801 7.0 0.400 0.25 1.80 0.050 51.0 189.0 0.99174 3.00 0.55 11.4 6
4045 6.8 0.120 0.31 5.20 0.045 29.0 120.0 0.99420 3.41 0.46 9.8 7
3651 7.9 0.460 0.40 10.10 0.168 19.0 184.0 0.99782 3.06 0.62 9.5 5
2649 6.8 0.260 0.43 11.75 0.045 53.0 198.0 0.99690 3.26 0.55 9.5 5
3930 7.5 0.570 0.02 2.60 0.077 11.0 35.0 0.99557 3.36 0.62 10.8 6
3287 11.9 0.380 0.51 2.00 0.121 7.0 20.0 0.99960 3.24 0.76 10.4 6
5123 7.1 0.250 0.28 1.20 0.040 31.0 111.0 0.99174 3.18 0.53 11.1 5
3195 8.3 0.615 0.22 2.60 0.087 6.0 19.0 0.99820 3.26 0.61 9.3 5

Al igual que sucede con los modelos clasificadores, hay multitud de modelos regresores, por lo que vamos a ir uno a uno para finalmente, poder escoger el mejor:

1.Regresión lineal múltiple (mínimos cuadrados ordinarios):

Al tratarse de una regresión lineal con multitud de variables, llamaremos al dato X.

Nuestro principal objetivo es minimizar el error cuadrático medio entre los resultados reales y los estimados.

In [89]:
from sklearn.linear_model import LinearRegression

linear_model = LinearRegression(fit_intercept=True, 
                                normalize=False,
                                n_jobs=1
                                )
In [90]:
# Hay que llamar al método .fit() para entrenar la regresión:

linear_model.fit(X=train_r["features"],
                y=train_r["label1", "quality"])
Out[90]:
LinearRegression(n_jobs=1)

Buscamos las Betas de cada una de las variables independientes junto con el intercepto, por lo que dispondríamos de una fórmula de regresión lineal:

In [91]:
print("beta 0 o intercept: " + str(linear_model.intercept_))
print("beta 1: " + str(linear_model.coef_[0]))
print("beta 2: " + str(linear_model.coef_[1]))
print("beta 3: " + str(linear_model.coef_[2]))
print("beta 4: " + str(linear_model.coef_[3]))
print("beta 5: " + str(linear_model.coef_[4]))
print("beta 6: " + str(linear_model.coef_[5]))
print("beta 7: " + str(linear_model.coef_[6]))
print("beta 8: " + str(linear_model.coef_[7]))
print("beta 9: " + str(linear_model.coef_[8]))
print("beta 10: " + str(linear_model.coef_[9]))
print("beta 11: " + str(linear_model.coef_[10]))
beta 0 o intercept: 55.21530587279736
beta 1: 0.06255635508077703
beta 2: -1.181690257998828
beta 3: 0.07173958235946717
beta 4: 0.034997598499811024
beta 5: -0.6622656077456244
beta 6: 0.00827382622474062
beta 7: -0.0026224020912718987
beta 8: -55.1369181500756
beta 9: 0.5999538056351338
beta 10: 0.8229176906381876
beta 11: 0.2738251492661459

Podemos medir el rendimiento de la regresión lineal simple tanto con el coeficiente de determinación (R cuadrado) como con el MSE o RMSE.

Tenemos el objetivo de saber cuanto nos hemos acercado a nuestro número real en la predicción (no en calcular cuántos acertamos), por lo que aplicaremos el método .score().

Coeficiente de determinación R cuadrado:

En esta ocasión lo haremos con el R cuadrado (coeficiente de determinación), un indicador que mide el porcentaje de la variabilidad del target explicado por la variabilidad de las features.

In [92]:
print("R cuadrado: " + str(linear_model.score(X=test_r["features"],
                                              y=test_r["label1", "quality"])))
R cuadrado: 0.27059302749465874

Hemos obtenido un R cuadrado muy bajo para nuestro conjunto de test, por lo que parece que no será el mejor tipo de métrica para nuestro problema de regresión. Ciertamente no es la medida de rendimiento más recomendada para las regresiones.

El R cuadrado es una medida de la relación lineal entre las features y el label a predecir, con un valor cercano a 1 indicaría que una larga proporción de la variabilidad de nuestra label está explicada o quitada gracias al modelo.

MSE (Error Cuadrádito Medio) y RMSE (Root Error Cuadrádito Medio):

Son otro tipo de métricas para medir el rendimiento de regresión, suelen ser las más utilizadas.

El MSE toma como primer argumento los valores reales del target, y como segundo los valores que queremos predecir.

In [93]:
from sklearn.metrics import mean_squared_error

# Aplicamos el MSE en el conjunto de testing
mse = mean_squared_error(y_true = test_r[["label1"]],
                         y_pred = linear_model.predict(test_r["features"])
                        )

# El RMSE únicamente es la raíz cuadrado del MSE.
rmse = np.sqrt(mse)

print("Nuestro MSE es: " + str(mse))
print("Nuestro RMSE es: " + str(rmse))
Nuestro MSE es: 0.5533326229331919
Nuestro RMSE es: 0.7438633093070204

El RMSE se puede entender en términos medios, fallamos las predicciones por 0.7438 de puntuación, sobre el valor real.

Parece que vamos a tener que aplicar regularizaciones. Estas consisten en alterar ligeramente la formulación matemática de un modelo de ML, con la intención de prevenir el sobreajuste.

Su principal finalidad, es añadir hiperparámetros a nuestro modelo que nosotros mismos elegimos, estos incluyen penalizaciones con la finalidad de que, cuando el modelo se entrene, lo haga con la intención de recudir el error.

El modelo ahora tendrá que conseguir un balance entre ajustar bien, y no estar demasiado penalizado.

Aplicaremos las regresión Ridge, pero tambien existen la regresión LASSO y Elastic Net:

2.Regresión Ridge:

Es la regresión lineal que aplica la regularización L2, que tienen la finalidad de que los parámetros estimados por nuestro modelo no tengan valores absolutos demasiado grandes.

In [94]:
from sklearn.linear_model import Ridge

# La formulamos:
ridge = Ridge(alpha= 4,  # Hiperparámetro de regularización, marca la intensidad de la contracción
                         # de los parámetros. Probamos con 4
             fit_intercept = True)

# Entrenamos:
ridge.fit(X=train_r["features"],
         y=train_r["label1", "quality"])
Out[94]:
Ridge(alpha=4)

Ahora hay que ver si los parámetros beta estimados se han hecho más pequeños, comparados con los de la regresión lineal simple por mínimos cuadrados anterior.

In [95]:
print("beta 0 o intercept: " + str(ridge.intercept_))
print("beta 1: " + str(ridge.coef_[0]))
print("beta 2: " + str(ridge.coef_[1]))
print("beta 3: " + str(ridge.coef_[2]))
print("beta 4: " + str(ridge.coef_[3]))
print("beta 5: " + str(ridge.coef_[4]))
print("beta 6: " + str(ridge.coef_[5]))
print("beta 7: " + str(ridge.coef_[6]))
print("beta 8: " + str(ridge.coef_[7]))
print("beta 9: " + str(ridge.coef_[8]))
print("beta 10: " + str(ridge.coef_[9]))
print("beta 11: " + str(ridge.coef_[10]))
beta 0 o intercept: 1.2381137230091213
beta 1: 0.004504318617346953
beta 2: -1.2645550090891147
beta 3: 0.07110076632488319
beta 4: 0.013812556331848714
beta 5: -0.4763463124564348
beta 6: 0.008256147742815289
beta 7: -0.0022887001355284082
beta 8: -0.03832715266051679
beta 9: 0.3106174313071383
beta 10: 0.6306129228393811
beta 11: 0.34220901758648564
In [96]:
print("Ridge: " + str(ridge.score(X=test_r["features"],
                                  y=test_r["label1", "quality"])))

# Aplicamos el MSE en el conjunto de testing
mse = mean_squared_error(y_true = test_r[["label1"]],
                         y_pred = ridge.predict(test_r["features"])
                        )

# El RMSE únicamente es la raíz cuadrado del MSE.
rmse = np.sqrt(mse)

print("Nuestro MSE es: " + str(mse))
print("Nuestro RMSE es: " + str(rmse))
Ridge: 0.2673275654482643
Nuestro MSE es: 0.5558098225588192
Nuestro RMSE es: 0.7455265404791563

Las betas estimadas se han visto reducidas y nuestro score se ha visto reducido a 0.2673 en el caso de la regresión Ridge, prácticamente el RMSE continúa igual, por lo que apenas hemos mejorado el modelo anterior.

De todas formas las regresiones lineales no suelen ser ni las más utilizadas ni las mejores, por lo que vamos a probar con otras.

3.Árboles regresores:

En este caso los árboles no son modelos lineales, por lo que no habrá que realizar transformaciones no lineales.

In [97]:
from sklearn.tree import DecisionTreeRegressor

# Creamos el GridSearch de árboles:
arbol_gs = GridSearchCV(estimator=DecisionTreeRegressor(),
                        param_grid={"max_depth": [2,3,5,7],
                                    "min_samples_split": [2,12,14,17],
                                    "min_samples_leaf": [2,12,14,17]},
                        scoring="neg_mean_squared_error", # MSE negativo
                        cv=5,
                        verbose=1
                       )

# Entrenamos, validamos y elegimos el mejor árbol:
arbol_gs.fit(X=train_r[["features"]],
             y=train_r["label1", "quality"]
            )
pass
Fitting 5 folds for each of 64 candidates, totalling 320 fits
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done 320 out of 320 | elapsed:    3.9s finished

Una vez ha impreso cómo ha ido entrenando y validando cada modelo con validación cruzada, podemos preguntarle cuál ha sido el mejor modelo en validación con el atributo:

In [98]:
arbol_gs.best_estimator_
Out[98]:
DecisionTreeRegressor(max_depth=5, min_samples_leaf=12)
In [99]:
arbol_gs.best_score_
Out[99]:
-0.5775095995289568

En este caso el MSE obtenido es negativo (lo hemos indicado en los hiperparámetros anteriormente). Es exactamente el mismo significado que si estubiese en positivo, solo que en vez de minimizar el MSE lo maximizamos.

Podemos extraer el MSE que hemos estado utilizando hasta ahora (positivo) para tener mayor facilidad en comparar.

In [100]:
# Medimos:
mse = mean_squared_error(y_true=test_r["label1", "quality"],
                         y_pred=arbol_gs.predict(test_r[["features"]]))
# Calculamos el RMSE:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
MSE: 0.5934165605402605
RMSE: 0.7703353558939512

Podemos realizar otro modelo de árboles probando de la siguiente manera, hasta aconseguir el mejor score:

In [101]:
from sklearn.tree import DecisionTreeRegressor

# Hay que hacer fórmula del árbol:
tree_reg = DecisionTreeRegressor(criterion="mse",
                                 max_depth=7,
                                 min_samples_split=8,
                                 min_samples_leaf=8,
                                )

# Entrenamos:
tree_reg.fit(X=train_r[["features"]],
             y=train_r["label1", "quality"]
            )

# Medimos:
mse = mean_squared_error(y_true=test_r["label1", "quality"],
                         y_pred=tree_reg.predict(test_r[["features"]]))
# Calculamos el RMSE:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
MSE: 0.6161051852180403
RMSE: 0.784923681142339

4. Ensembles:

Tambien existen ensembles de regresores, tal como hemos podido verlos en las clasificaciones.

  • En el caso de los ensembles de Bagging la técnica es hacer la media de las predicciones de cada modelo lanzado.

  • En el caso de los ensembles de Boosting, se va corrigiendo el MSE de los modelos anteriores con los nuevos.

Vamos a hacer un Random Forest Regresor:

In [102]:
from sklearn.ensemble import RandomForestRegressor

rf_reg = RandomForestRegressor(n_estimators=40, # 40 árboles
                               criterion="mse",
                               max_depth=7, 
                               min_samples_split= 12,
                               min_samples_leaf=8,
                               bootstrap=True
                              )

# Entrenamos:
rf_reg.fit(X=train_r[["features"]],
           y=train_r["label1", "quality"]
          )

# Medimos:
mse = mean_squared_error(y_true=test_r["label1", "quality"],
                         y_pred=rf_reg.predict(test_r[["features"]]))
# Calculamos el RMSE:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " + str(rmse)) 
MSE: 0.5089146038166944
RMSE: 0.7133825087683987

O bien podemos hacerlo vía GridSearchCV para seleccionar los mejores hiper parámetros:

In [103]:
rf_reg_2 = GridSearchCV(estimator=RandomForestRegressor(),
                        param_grid={"n_estimators": [150],
                                    "max_depth": [2,12,14,17,25],
                                    "max_features": ["sqrt", 3, 4]},
                        scoring="neg_mean_squared_error", # MSE negativo
                        cv=5,
                        verbose=1,
                        n_jobs=-1
                       )
# Entrenamos, validamos y elegimos el mejor árbol:
rf_reg_2.fit(X=train_r[["features"]],
             y=train_r["label1", "quality"]
            )
pass
Fitting 5 folds for each of 15 candidates, totalling 75 fits
[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed:   15.6s
[Parallel(n_jobs=-1)]: Done  75 out of  75 | elapsed:   29.0s finished
In [104]:
rf_reg_2.best_estimator_
Out[104]:
RandomForestRegressor(max_depth=17, max_features=4, n_estimators=150)
In [105]:
rf_reg_2.best_score_
Out[105]:
-0.4740435943496492
In [106]:
# Medimos:
mse = mean_squared_error(y_true=test_r["label1", "quality"],
                         y_pred=rf_reg_2.predict(test_r[["features"]]))
# Calculamos el RMSE:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " + str(rmse))
MSE: 0.4726078258566267
RMSE: 0.6874647815391176

5.Regresión por Nearest Neighbors:

Es un modelo no lineal, el concepto es como el del clasificador, mira los vecinos cercanos a los puntos a predecir, y hace la media de los valores de dichos vecinos.

In [107]:
from sklearn.neighbors import KNeighborsRegressor

knn_reg = KNeighborsRegressor(n_neighbors=7,
                              metric="euclidean")

# Entrenamos:
knn_reg.fit(X=train_r[["features"]],
            y=train_r["label1", "quality"]
           )

# Medimos:
mse = mean_squared_error(y_true=test_r["label1", "quality"],
                         y_pred=knn_reg.predict(test_r[["features"]]))
# Calculamos el RMSE:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " + str(rmse)) 
MSE: 0.661865121988645
RMSE: 0.8135509338625609

O bien podemos hacerlo vía GridSearchCV para seleccionar los mejores hiper parámetros:

In [108]:
knn_gs = GridSearchCV(estimator=KNeighborsRegressor(),
                      param_grid={"n_neighbors": [1,2,3,5,9,13,45]},
                      scoring="neg_mean_squared_error",
                      cv=5)

# Entrenamos, validamos y elegimos el mejor árbol:
knn_gs.fit(X=train_r[["features"]],
             y=train_r["label1", "quality"]
            )
pass
In [109]:
knn_gs.best_estimator_
Out[109]:
KNeighborsRegressor(n_neighbors=13)
In [110]:
knn_gs.best_score_
Out[110]:
-0.6757984450786749
In [111]:
# Medimos:
mse = mean_squared_error(y_true=test_r["label1", "quality"],
                         y_pred=knn_gs.predict(test_r[["features"]]))
# Calculamos el RMSE:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " + str(rmse)) 
MSE: 0.6435189304622504
RMSE: 0.8021963166596132

6.Suport Vector Regression (SVR):

Al igual que su clasificador (SVM), el cual es muy conocido para clasificación, tambien existen en regresión, pero no son tan destacables como en clasificación:

In [112]:
from sklearn.svm import SVR

svr = SVR(C=0.5,
          kernel="poly", # Kernel polinómica de tercer grado
          degree=3
         )

# Entrenamos
svr.fit(X=train_r[["features"]],
        y=train_r["label1", "quality"]
       )

# Medimos:
mse = mean_squared_error(y_true=test_r["label1", "quality"],
                         y_pred=svr.predict(test_r[["features"]]))
# Calculamos el RMSE:
rmse = np.sqrt(mse)

print("MSE: " + str(mse))
print("RMSE: " + str(rmse)) 
MSE: 0.7166403063905004
RMSE: 0.8465461041139463

Hemos entrenado unos cuantos modelos, por lo que a partir de estos, seleccionaremos cual ha sido el mejor de todos para nuestra predicción regresora sobre la calidad del vino.

In [113]:
# Metemos los modelos en los que hemos utilizado el GridSearchCV en una lista
modelos_entrenados = [arbol_gs, rf_reg_2, knn_gs]

# Generamos una lista para meter los .best_score_ y los metemos:
lista_scores = []

for modelo in modelos_entrenados:
    lista_scores.append(modelo.best_score_)
In [114]:
lista_scores
Out[114]:
[-0.5775095995289568, -0.4740435943496492, -0.6757984450786749]

Ya tenemos un modelo ganador, en este caso el modelo seleccionado es el rf_reg_2, es decir el Random Forest.

Nuestro modelo vencedor tiene un RMSE: 0.687322064531013, por lo que las predicciones del modelo final se alejan en promedio 0.6873 unidades del valor real, bastante poco. Es el RMSE más bajo de todos los modelos entrenados, por lo que es el claro vencedor al disponer del menor error cuadrático medio.

El modelo resultante es el siguiente:

In [115]:
rf_reg_2.best_estimator_
Out[115]:
RandomForestRegressor(max_depth=17, max_features=4, n_estimators=150)

De este modo, podemos dar por concluída la búsqueda del mejor modelo regresor para predecir la calidad del vino a partir de las features facilitadas en nuestro dataset.